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ABSTRACT 

We report results from numerical simulations of star formation in the early universe that focus 
on the dynamical behavior of metal-free gas under different initial and environmental conditions. 
In particular we investigate the role of turbulence, which is thought to ubiquitously accompany the 
collapse of high-rcdshift halos. We distinguish between two main cases: the birth of Population III. 1 
stars - those which form in the pristine halos unaffected by prior star formation - and the formation 
of Population III. 2 stars - those forming in halos where the gas has an increased ionization fraction. 
We find that turbulent primordial gas is highly susceptible to fragmentation in both cases, even for 
turbulence in the subsonic regime, i.e. for rms velocity dispersions as low as 20 % of the sound speed. 
Fragmentation is more vigorous and more widespread in pristine halos compared to pre-ionized ones. 
If such levels of turbulent motions were indeed present in star- forming minihalos, Pop III. 1 stars 
would be on average of somewhat lower mass, and form in larger groups, than Pop III. 2 stars. We 
find that fragment masses cover over two orders of magnitude, suggesting that the Population III 
initial mass function may have been much broader than previously thought. This prompts the need 
for a large, high-resolution study of the formation of dark matter minihalos that is capable of resolving 
the turbulent flows in the gas at the moment when the baryons become self-gravitating. This would 
help to determine the applicability of our results to primordial star formation 



Subject headings: stars: formation 
of state 



1. INTRODUCTION 



stars: mass function - early universe - hydrodynamics - equation 



Over the course of the last decade, work by a number 
of groups has led to the development of a widely ac- 
cepted picture for the formation of the first stars, the so- 
called Population III (or Pop. Ill) stars. In this picture, 
the ve ry first stars (Popul atio n III. 1 in the nom encla- 
ture of |Tan fe McKee|2008| and |O'Shea et~aTp 008) form 
within small dark matter halos that have total masses 
M ~ 10 6 M , virial temperatures of around a thou- 
sand K, a nd that are assembled at redshifts z ~ 25-30 
or above (|Bromm fe Larson 2004 Glover 2005 Bromm 



et al.||2005 r Gas falling into one of these small dark 
matter halos is shock-heated to a temperature close to 
the virial temperature of the halo, and thereafter cools 
via H2 rotational and vibrational line emission. The H2 
that enables the gas to cool is primarily formed via the 
gas-phase reactions 

H + c-^H-+7, (1) 

H-+H^H 2 +e", (2) 

where the required free electrons are those that remain 
in the gas after cosmological recombination at z ~ 1100. 
The low abundance of these electrons, plus the fact that 
the gas starts to recombine further once it collapses into 
the dark matter halo, limits the amount of H2 that can 
form in this manner; the typical fractional abundances 
of H 2 that result are a few times 10 -3 . The limited H 2 
abundance, the low critical density of the H2 molecule, 



and the large energy separation of its lowest accessible 
rotational levels combine to significantly limit the ex- 
tent to which the gas can cool. The minimum tempera- 
ture reached depends upon the dynamical history of the 
collapse, but is typically around 200 K. This minimum 
temperature is reached at a density of around 10 4 cm~ 3 , 
comparable to the critical density of H 2 (the density at 
which its level populations reach their local thermody- 
namical equilibrium (LTE) values), and at higher densi- 
ties, the gas begins to reheat. Gas falling into the dark 
matter halo therefore tends to fragment once it reaches 
this temperature and density, and the resulting frag- 
ments have masses of the order of the local Jeans mass, 
Mj ~ 1000 Mq. 1 Following this period of fragmentation, 
the gas in the fragments undergoes a further period of 
gravitational collapse, eventually reaching a density of 
n ~ 10 8 -10 10 cm -3 , at which point three-body H 2 for- 
mation becomes effective. This rapidly converts almost 
all of the available atomic hydrogen into H 2 , but the 
significant energy release that accompanies this process 
heats the gas, which typically attains a temperature of 
1000-2000 K during this phase of the collapse. 
Continued collapse next leads to the gas becoming 

Note that the term 'fragmentation' is often used in the liter- 
ature even in the case where only one object forms at the center 
of the clump, even though, strictly speaking, the term refers to a 
system separating into several distinct parts with separate evolu- 
tionary paths. 
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optically t hick in the main H? cooling lines (at 
10 10 cm" " 



Ripamonti fc Abel|2004[|Yoshida et al.|2006| ) 



the onset o f collision- induced emissi on cooling (at n 
10 14 cn 



Ripamonti & Abel 2004), and finally to the 



gas becoming optically thick in the continuum as well 



as in the H 2 lines (n ~ 10 cm" , Yoshida, Omukai 
& Hernquist 2008). H 2 dissociation cooling allows for 
further collapse, but once all of the H2 is gone, the col- 
lapse becomes fully adiabatic. Detailed simulations of 
the collapse, using adaptive mesh refinement (AMR; see 
e.g. Abel, Bryan, fc Norman||2002| |Q'Shea fc Norman 
2007D or smoothed particle hydrodynamics (SPH; see 
e.g. IBromrri fc Loeb 20041 1 Yoshida et al.||2006| |Yoshida,| 
Omukai fc Hernquist|2008 ) show little or no evidence for 



ragmentation between n ~ 10 4 cm -3 and the onset of 
this adiabatic evolution, and it is natural to associate 
the latter with the formation of a primordial protostar. 
The amount of gas incorporated into the protostar at 
this point is small, approximately 0.01 M , but it is sur- 
rounded by a massive, dense envelope and hence begins 
to accrete rapidly. If the envelope does not fragment, and 
if accretion onto the protostar is unhindered by radiative 
feedback, then the final mass of the star can be very 
large. For example, O'Shea & Norman (2007) estimate 
that its mass could lie anywhere in the range M ~ 20- 
2000 M Q , depending on the environment and dynamical 
history of the gas. The degree to which feedback will sup- 
press accretion remains uncertain, but the most effective 
potential feedba ck mechanisms, Lym an-q scattering and 
photoionization ( McKee & Tan 2008 ) become important 
only once the mass of the star exceeds 20 M . The ex- 
pectation is therefore that all Population III. 1 stars were 
massive, with masses typically of the order of 100 M Q 
or more, thereby also explaining why none of these stars 
appear to have survived until the present day. If a Pop- 
ulation III. 1 star of this mass does form, then it will 
rapidly dissociate H 2 throughout the halo, thereby sup 



pressing further star formation (Omukai & Nishi 1999 - 
Glover fc Brand|2001 ). This fact, together with the lack 



of fragmentation seen in the simulations, is often taken 
to imply that only a single Population III. 1 star will form 
in each dark matter halo. 

This widely accepted picture for Population III star 
formation also provides for a second mode of Pop. Ill 
star formation. This occurs in metal-free gas which has 
been ionized by radiation from a previous generation of 
Population III stars. Following the death of these stars, 
the gas recombines, and the elevated fractional ioniza- 
tion in the recombining gas allows more H2 to form. The 
enhanced H2 fraction enables the gas to cool to a lower 
temperature, which in turn increases the effectiveness of 
cooling by the singly-deuterated hydrogen molecule, HD. 
In the standard Pop. III. 1 scenario, HD cooling is of lim- 
ited importance, but in this sec ond scenario, termed P op- 
ulation III.2 star formation by |Tan fc McKee| p008| , it 
becomes dominant, cooling the gas down to the temper- 
ature of the cosmic microwave backgro und (|Nagakura fc| 
Omukai 2005 ; Johnson & Bromm 2006 ; Yoshida, Omukai 



fc Hernquist|2007 ). The lower gas temperature, plus the 
higher critical density of HD (n cr it ~ 10 6 cm~ 3 , compared 
with n cr it ~ 10 4 cm~ 3 for H2), means that fragmentation 
occurs somewhat later, and produces significantly lower 
mass fragments, with characteristic masses M ~ 1OOM . 
Subsequently, the evolution of these fragments is believed 



to proceed in the same fashion as described above, result- 
ing ultimately in the formation of a primordial protostar 
of similar mass. However, since this protostar is embed- 
ded in a less massive envelope, with a lower accretion 
rate, the final stellar mass is thought to be an order of 
magnitude or so smaller, M ~ 10 M Q . Nevertheless, 
this still corresponds to what we would call, by Galac- 
tic standards, a massive star, and moreover one which is 
more than capable of dissociating H 2 throughout a large 
volume of the halo, thereby suppressing further star for- 
mation. 

In the past few years, it has become possible to par- 
tially test this picture using numerical simulations that 
begin with the proper cosmological initial conditions and 
follow the collapse of the gas all the way to protostellar 
densities. These simulations confirm that fragmentation 
during the initial collapse phase is ineffective: in gen- 
eral, only a single protostar forms, although in roughly 
20% of cases, the gas fragments into two clum ps, form- 
ing a wide binary ( |Turk, Abel fc O'Shea|[2009l M. Turk, 
private communication) . However, most of these stud- 
ies were prevented by technical limitations from follow- 
ing the evolution of the infallmg gas after the formation 
of the first protostar. The technical problem involves 
the hydrodynamical timestep: as the collapse is followed 
down to protostellar densities, this becomes extremely 
short (At ~ 10~ 3 yr), making it computationally infeasi- 
ble to follow the evolution of the surrounding dense gas 
over any significant time period. 

In studies of contemporary star formation, a s imilar 
problem is avoided by the use of sink particles ( |Bate,| 
Bonnell & Price 1995[ ). When gravitationally bound re- 
gions of gas collapse below the scale on which they can 
be spatially resolved, they are replaced in the simulation 
by a sink particle of the same mass. This can accrete gas 
from its surroundings, but otherwise interacts only via 
gravity. Recent simulations of Population III star forma- 
tion that have used sink particles to allow the evolution 
of the gas to be followed beyond the formation of the 
first protostar consistently find evidence for more exten- 
sive fragmentation than is assumed in the conventional 
picture of Population III star formatio n ( Clark, Glover fc 
Klessen||2008f |Stacy, Greif fc Bromm||20lo| Clark etaT, 
in prep.). If fragmentation is as common as these results 
suggest, and most Population III stars form in multiple 
systems rather than as single stars, then this has pro- 
found implications for the final masses of the stars, their 
production of ionizing photons and metals, the rate of 
high redshift gamma ray bursts, and many other issues. 

It is therefore important to better understand the phys- 
ical basis of fragmentation in these systems. To do this, 
we would ideally like to have a large sample of sim- 
ulations, as basing our arguments on only one or two 
realizations of Population III star formation leaves one 
open to the possibility that these realizations may not 
be typical. Unfortunately, simulations starting from cos- 
mological initial conditions and following collapse all the 
way to protostellar densities are computationally costly, 
making it difficult to explore a large region of parameter 
space. In these simulations, it is also difficult to be cer- 
tain about which aspects of the included physics are the 
most important for driving fragmentation. Simulations 
that start from simpler initial conditions and that focus 
on exploring the importance of a single free parameter 
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therefore play an important role, which is complemen- 
tary to that of fully cosmological simulations. A good 
exam ple is the recent work by Machida a nd collabora- 



tors ( |Machida||2008| |Machida et al.||2009 ) which exam- 
ined the influence of the initial rotational energy of the 
gas, and found that zero metallicity clouds with suffi- 
cient initial rotational energy could fragment into tight 
binaries. Another example, and one which explores a 
mu ch lower density regime, is the study by |Jappsen et| 



al. (2009ajb), who showed that the fragmentation behav- 
ior depends sensitively on the adopted density profile of 
the primordial halo. 

In this paper, we perform a similar study into the ef- 
fects of the initial turbulent energy of the gas. Our mo- 
tivation comes from the fact that recent high-resolution 
cosmological simulations have shown that the self- 
gravitating regions in which Pop. Ill star formati on oc- 
curs contain significant turbulent motions (see e.g. Greif 



et al. 2008). From studies of present-day star forma- 
tion, we know that such turbulent motions can lead to 
the fra gmentation of initially marginally unstable gas 



clouds (Klesscn 2001 Goodwin et al. 2004a b Mac Low 



& Klesscn 2004[ |Attwood et al.|2009|) , and so it is mstruc 



tive to investigate how such conditions could affect the 
'standard' picture of primordial star formation. In a Pop. 
Ill analogue of the studies on present-day star formation, 
we will look at how subsonic turbulence affects the col- 
lapse of marginally super-critical Bonnor-Ebert spheres, 
in an attempt to quantify the effects of turbulently in- 
duced fragmentation on the mass function of Pop. Ill 
stars. Central to this study will be the use of the sink 
particles to allow us to capture the evolution of the gas 
cloud beyond the collapse of the first region. 

The paper is laid out in the following manner. Our 
modifications t o the cosmological SPH code Gadget 2 
( Springel|20 05) are outlined in f|2](with further details of 
the chemical networks and heating and cooling described 
in the Appendix). The initial conditions for this numer- 
ical experiment are described in £j3| including the defini- 
tions of Pop. III. 1 and III. 2 that are used in our study. 
The details of the fragmentation process seen in the sim- 
ulations are discussed in { 4] and the long-term chemical 
and thermodynamical evolution of the infalling envelope 
are described in S|5j We discuss the implications of this 
study for the initial mass function (IMF) of primordial 
stars in Sj6j and summarize the main points of this paper 

in CD 

2. NUMERICAL METHOD 

We model the evolution of the gas in our simulations 
using a modified version of the Gadget 2 smoothed parti- 
cle hydrodynamics code (Springcl 2005). We have modi- 
fied the publicly available code in several respects. First, 
we have added a s ink particle implementation, based on 
the prescription in Bate, Bonnell & Price ( 1995 ), to allow 
us to follow the evolution of the gas beyond the point 
at which the first protostar forms. Our particular im- 
plementation is deriv ed from the one first described in 
Jappsen et al. (2005). We briefly describe the ideas be- 
hind the algorithm here. The actual numerical values 
used for the parameters discussed here are given later, in 

Sink particles are not really added to the code in the 
sense that a new particle is introduced, but instead a 



normal SPH particle is turned into a sink particle once 
certain criteria have been met. The particle undergoes 
a series of tests once it has reached a threshold density. 
The first is to see whether the candidate SPH particle 
is sufficiently far away from any other sinks, measured 
in terms of the sink particle's accretion radius, r acc . We 
adopt a conservative value of 2r acc . Next, we check to 
see whether the smoothing length of the particle is less 
than the 'accretion radius' of the sink particle that it 
will become. This ensures that when the sink particle 
forms it can instantly accrete at least ~ 50 neighboring 
particles (we adopt 50 neighbors in these simulations). 
The third test is to make sure that the candidate sink 
particle and its neighbors are on the same integration 
time-step. 

Once these three preliminary criteria are met, the dy- 
namical state of the possible sink particle and its neigh- 
bors are assessed to ensure that the particles are indeed 
undergoing gravitational collapse and are not about to 
re-expand from their dense state. This takes the form of 
a further four tests. Firstly, we require that 



a < 



1 



(3) 



where a is the ratio of the thermal energy to the mag- 
nitude of the gravitational energy of the particles. Sec- 
ondly, we ensure that 



a + (3 < 1 



(4) 



where j3 is the ratio of rotational energy to the magni- 
tude of the gravitational energy. The third condition is 
that the total energy of the particles must be negative 
(which actually renders the above checks redundant, but 
can help to improve the computational efficiency of the 
code) . Finally, the forth test is that the divergence of the 
accelerations must be less than zero. This final check en- 
sures that the group of particles is not in the process of 
being tidally disrupted or bouncing. If all of these tests 
are passed, the particle becomes a sink and the mass and 
linear momentum of the neighbors are added to the sink 
particle. 

As the simulation progresses, the sink particles are 
then allowed to accrete other gas particles that fall within 
the accretion radius. As in the Bate et al. (1995) pre- 
scription, several tests must be passed before any SPH 
particles can be accreted by the sink. First, it must ob- 
viously be bound and moving toward the candidate sink. 
Second, in the case where there is more than one sink 
present, it must be more bound to the candidate sink 
than any other sink in the simulation. Finally, the SPH 
and sink particles need to be on the same integration 
time-step (to ensure temporal momentum conservation) . 
Once these conditions are passed, the mass and linear 
momentum of the SPH particle are added to the sink 
particle. Note that accretion onto existing sink particles 
is done before any new candidate sinks are considered. 

In addition to the sink particles, we have also imple- 
mented an external pressure term (e.g. Benz||1990 ), that 
enables us to model a constant pressure boundary, as 
opposed to the vacuum or periodic boundary conditions 
that are the only choices available in the standard ver- 
sion of Gadget. We modify the standard gas pressure 
contribution to the Gadget 2 momentum equation, 



4 



dt 



3 



pi 



Pj 



by replacing p and Pj with Pj — P cxt and Pj 



(5) 

Pcxt 



respectively, where P ox t is the external pressure, and 
all quantities have the u sual m eaning, consistent with 
those used by Springel (2005). The pair- wise nature 
of the force summation over the SPH neighbors ensures 
that P ext cancels for particles that are surrounded by 
other particles. At the edge, where the term does not 
disappear, it mimics the pressure contribution from a 
surrounding medium. The values of P ex t used in this 
study are 3 x 10 7 k^K cm -3 for the Pop III. 1 clouds, and 
7.5 x 10 6 fce K cm~ 3 for the Pop III. 2 clouds, where k^ is 
the Boltzmann constant. These confining pressures are 
similar to the internal pressure of the clouds in each case 
(see the temperatures and densities in q3] below). 

Finally, we have added a treatment or primordial gas 
chemistry and cooling to the code. To model the thermal 
evolution of the gas, we use an operator-split formalism, 
which treats the effects of radiative and chemical heating 
and cooling separately from compressional heating. The 
influence of radiative and chemical heating and cooling 
on the thermal energy of each particle can thus be formu- 
lated as an ordinary differential equation (ODE), which 
can then be solved simultaneously with the chemical rate 
equations (also ODEs) that describe the chem ical evolu- 
tion of the gas. As |Turk, Abel fc O'Shea| ( |2009[ ) have pre- 
viously discussed, the strong coupling between the chem- 
ical and thermal evolution of high density primordial gas 
that results from the importance of three-body H2 for- 
mation heating and H2 collisional dissociation cooling 
renders it vital to treat cooling and chemistry simultane- 
ously; treatments that do not do so will produce numeri- 
cally stable results only with great difficulty. Full details 
of our chemical model and cooling function are given in 
the Appendix. 

3. INITIAL CONDITIONS 

Since the aim of this study is to investigate the frag- 
mentation properties of primordial gas, we use initial 
conditions that allow us to perform a controlled numer- 
ical experiment. The clouds start as unstable Bonnor- 
Ebert (BE) spheres into which we inject a subsonic tur- 
bulent velocity field. The BE sphere is made by allowing 
a gravitationally stable cloud of gas to evolve under its 
own self-gravity until the system has settled into a sta- 
ble, centrally condensed configuration. During this ini- 
tial settling phase, the gas temperature is held constant 
and the chemical evolution is not followed. On the scales 
of interest in this study, gravitational forces from the 
dark matter are negligible compared to the self-gravity 
of the gas, and so for simplicity, we do not include dark 
matter in our models. 

By re-scaling the mass and temperature of the cloud, 
we are then able to choose initial conditions that are 
gravitationally unstable, and which are similar to the ini- 
tial conditions for the Pop. III. I and Pop. III. 2 star for- 
mation channels that occur in cosmological simulations. 
All of the simulations presented in this paper start with 
a maximum central cloud number density of 10 5 cm~ 3 . 
The Pop. III. I simulations start with an initial temper- 



ature of 300K and contain 1000 M of gas. For the Pop. 
III. 2 simulations, we adopt an initial temperature of 75 K 
and examine two different initial masses. The first set 
of Pop. III. 2 simulations contain 150 M Q of gas, chosen 
such that they have the same ratio of thermal to grav- 
itational energy as the Pop. III. 1 simulations, roughly 
0.15. As such, these Pop. III. 2 simulations contain the 
same initial number of Jeans masses of gas as the Pop. 
III. I simulations. For a uniform sphere, the number of 
Jeans masses is given by (Ptherm/|P g rav|)~ 3 / 2 , so these 
clouds have roughly 3 Jeans masses in the initial con- 
figuration. All else being equal, if these clouds were to 
evolve isothermally from this point on, they would have 
the same propensity to fragment as the Pop. III. 1 clouds, 
since this is directly related to the initial ratio of grav- 
itational to thermal energy. Note that we have chosen 
a slightly subvirial configuration for our initial setup, to 
ensure that the clouds are still able to collapse when 
the turbulent motions are included. In addition to these 
simulations, we also performed a second set of Pop. III. 2 
simulations that start with the same gas mass as in the 
Pop. III.l case (1000 M Q ). In this case, the Pop. III.2 
simulations are initially more Jeans unstable than their 
Pop. III.l counterparts, and therefore might be expected 
to fragment significantly more. 

In our Pop. III.l simulations, we set the initial frac- 
tional abundances of H2, H + , HD and D + to ih 2 = 10~ 3 , 
aiH+ = 10~ 7 , ihd = 3 x 10~ 7 and x B + = 2.6 x 10~ 12 , re- 
spectively. Our values for ih 2 , Xh+ an( i ^hd are typical 
of the values found at these densities in cosmo logical sim- 
ulatio ns of Pop. III.l star formation (see e.g. |Greif et al 



2008 1 , and account for the fact that the HD /H2 ratio is 
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In the case of D + , fractionation 
is unimportant at our starting temperature, and so we 
simply set D+/H+ = 2.6 x 10 -5 . In our Pop. III. 2 simula- 
tions, we adopt the same initial H + and D + abundances, 
but set xh ? = 3 x 10 -3 and xhd = 3 x 10 -6 , following 
Greif et al.| ( |2008[ ). In both the Pop. III.l and Pop. 111.2 
simulations, we assumed that all of the helium remained 
neutral, and set the initial abundances of all of our other 
tracked species to zero. 

Within the BE spheres, we impose a turbulent veloc- 
ity field that has a power spectrum of P(k) oc k~ 4 . We 
assume that the turbulence considered here has its ori- 
gin in gravitationally driven flows that arise as the gas 
and dark matter virialize in mini-halos (Wise & Abel 
2007; Greif et al. 2008; Klessen & Hennebelle 2010). 
As the gas is compressible in nature the turbulent ve- 
locity field will have a power spectrum that is somewhat 
steeper than the standard Kolmogorov (1941) descrip- 
tion for incompressible flows. However, this deviation 
is small and we note that the the ability of a cloud to 
fragment is only weakly dependent on the power spec - 
trum of the turbulence (Delgado-Donate et al. 2004[ ). 
The three-dimensional root-mean-squared velocity in the 
turbulent field - which we will refer to as Aw tU rb ~~ is then 
scaled to some fraction of the sound speed c s in the ini- 
tial conditions. For the simulations presented here we 
use four different rms velocities: 0.1, 0.2, 0.4 and 0.8 c s . 
For an isothermal sound speed and an adiabatic index 
of 7 = 5/3, the corresponding ratios of the turbulent 
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Fig. 1. — Column density images showing the state of the Pop. III. 1 clouds after they have converted 10 percent of their mass (100 Mq) 
into sinks. The sink particles are denoted in the images by the white dots, and in each case we centre the image on the first sink particle to 
form in the simulation. For the run with the lowest level of turbulence (Ai) tur i, = 0.1c a ) the initial turbulent velocity field provides the gas 
with enough angular momentum to produce a small disk around the sink particle, but does not cause fragmentation of the gas. However, 
once the strength of the initial turbulent velocities is increased to as little as 0.2 times the initial sound speed in the gas, the turbulence 
induces fragmentation. At the point at which the simulations are shown here, the 0.1, 0.2, 0.4 and 0.8 c s runs have produced 1, 5, 31 and 
15 sink particles, respectively. Note that the sink particles in this study have an accretion radius of 20 au. 



to thermal energy are given by l/3(Avturb/c s ) 2 , yield- 
ing 0.0033, 0.0133, 0.0533 and 0.2133 respectively, for 
our chosen values of c s . In order to focus on the effects 
of the turbulence, we do not include any ordered rota- 
tion of the initial gas cloud. Note, however, that this 
does not imply that the initial angular momentum of 
the cloud is zero, since the imposed turbulent velocity 
field gives the cloud a small amount of angular momen- 
tum. Note that we only consider subsonic turbulence in 
this study since our clouds have only a few Jeans masses, 
and supersonic turbulence would unbind them. To study 
the effects of supersonic turbulence, one would have to 
look at clouds that are initially more Jeans unstable than 
those we study here. 
The clouds in our study are all modeled using 2,000,000 



SPH particles. Although this means that the mass res- 
olution is higher in the 150 M Pop. III. 2 clouds than 
in the other simulations, the Jeans mass at the point 
where sink particles form is well resolved in every case 
(see below). For this study, it is more important that the 
turbulence in all simulations is evolved with the same res- 
olution, hence our choice of a constant particle number 
throughout. 

Sink particles are created once the number density of 
the gas reaches 10 13 cm -3 , at which point the gas has a 
temperature of around 1200K. The corresponding Jeans 
mass at this density and temperature is 0.08 M Q . Our 
mass resolution in the Pop. III. 1 and the 1000 M Pop. 
III. 2 clouds is 27V nci gh77ip ar t = 0.05 M , where iV ncig h is 
the number of neighbors employed for force evaluations 
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Fig. 2. — The top panel shows the mass functions from those 
Pop. III. 1 simulations in which fragmentation occurs. In all cases 
the mass function is plotted at the point at which the total mass 
of gas converted to sink particles is IOOMq. Note that as accretion 
is ongoing, and the system is still young (t ~ 1000 yr), these will 
often not be the final masses of the sinks. The mass functions in the 
individual simulations differ substantially, although the combined 
mass function, shown in the bottom panel, exhibits a broad and 
flat distribution between masses of 0.4 and 20 Mq. 



(in our case 50), and TO par t is the mass of an SPH par- 
ticle. The 150 M Pop. III. 2 simulations have a mass 
resolution of 0.008 M Q . Once the candidate particle has 
passed the criteria described in Bate et al. (1995), it is 
replaced by a sink particle that can accrete gas parti- 
cles that fall within its accretion radius r acc , which we 
fix at 20 AU. Note that this radius is significantly larger 
than the Jeans radius at the point that the gas reaches 
the density threshold for sink creation, which is around 
6 AU. As such, all the fragmentation that we capture in 
this study is well resolved and not lying close to the lim- 
its of our resolution; if any 'artificial' fragmentation were 
to occur, the sink particle would immediately swallow 
the offending region and replace it with a single accret- 
ing point. We also prevent sink particles from forming 
within 2 r acc of one another. This prevents the forma- 
tion of sinks out of gas that in reality would by accreted 
by a neighboring sink particle before it could go into di- 
rect collapse by itself. Lastly, gravitational interactions 
between sinks, and between the sinks and the gas, are 



softened in the standard way in Gadget 2, using a fixed 
softening parameter of 5 AU for the sinks, and a variable 
softening parameter for the gas particles that is equiva- 
lent to their smoothing length. 

We note that the sink algorithm employed here does 
not allow sink particles to merge, in contrast to the 
'sticky' sinks introduced in Bromm et al. (2002). As 
such, our simulations may be biased towards low masses. 
Small, secondary sinks could be driven towards more 
massive ones via dynamical friction and coalesce with 
them. This effect, which acts to reduce the resulting 
number of fragments, has been seen in AMR simula- 
tions of massive star forming regions in the present-day 
universe (Krumholz et al. 2007, 2009). However, given 
the size of our sink particles, it is unclear whether sink- 
merging is really applicable: the protostellar radius is at 
most expe cted to be 0.5 au (and then only for a short 
time, e.g. Hosokawa & Omukai 2009), which is signif- 
icantly smaller than both our adopted accretion radius 
and the gravitational softening. In fact, in the simula- 
tions from present-day star formation in which extremely 
small sink particles are used - along with little, or no soft- 
ening—emerging events are found to be fairly rare (e.g. 



Bate 



2009) 



They do however occur, and this should be 
borne in mind when interpreting the results of this study. 
Also, given the size of our sink particles, we are unable 
to resolve any of the disks that would invariably form 
around the young protostars and it has been suggested 
that such structures would be unstable to fragmentation 
on these scales (e.g. 



Machida et al. 



r p008| ). Finally, our 
calculations do not include any model for the feedback 



proc esses that accompan y the formation of a young star 
(e.g. McKee k. Tan|2008 ) . Given these uncertainties, the 
sink particles are unable to say exactly what the shape 
of the final IMF will be, but rather measure how the gas 
can fragment at a given scale (our resolution), and how 
these fragments are likely to evolve, assuming feedback 
processes play only a minor role over the timescales in- 
vestigated in this study. These caveats should be borne 
in mind when interpreting the results in the following 
sections. 

4. THE FRAGMENTATION OF PRIMORDIAL GAS 

In our comparison of the different simulations per- 
formed in this study, we will compare the properties of 
the clouds after roughly 10 percent of their mass has 
been accreted (although in some cases we will mention 
in passing what happens as the simulations are advanced 
further) . In terms of looking at the ability of the clouds 
to fragment, comparing the different simulations at this 
point in their evolution ensures that they have turned 
the same faction of their initial number of Jeans masses 
into the sink particles. The exception in our analysis is 
the comparison between the 1000 M Q Pop. III. 1 and III. 2 
simulations, since in this case the Pop. III. 2 clouds are 
initially more Jeans unstable. We now go on to describe 
the fragmentation in the Pop. III. 1 and III. 2 channels in 
some detail. 

4.1. Pop. 111.1 clouds 

The panels in Fig. [I] show the column density distribu- 
tion in the Pop. III. 1 simulations after 10 percent of their 
mass has been accreted onto the sink particles. The im- 
ages show the inner 1300 AU of the cloud, centred on the 
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Fig. 3. — In the left-hand panel we show the mass evolution of the sink particles that form during the first 1200 years in the Pop. 
III. 1 simulations, during which just over 10 percent of the initial gas mass is accreted. The solid (black) lines chart the mass of individual 
sink particles while the dashed (blue) lines show the mass evolution of the entire cluster of sink particles. The right-hand panel shows the 
associated mass accretion rate of the cluster, both as measured by summing over all sink particles (blue dashed) and from an estimate 
based on the radial mass infall profile just before the formation of the first sink particle (red dashed). 



first sink particle to form in each simulation. We can see 
that the clouds with Av tU rb > 0.2c s fragment into small 
clusters of sink particles; the runs with Ai> tur b = 0.2c s , 
0.4c s and 0.8c s form 5, 31 and 15 sink particles respec- 
tively. These calculations clearly demonstrate that tur- 
bulent subsonic motions are able to promote fragmen- 
tation in primordial gas clouds. Only in the case with 
A^turb = 0.1c s does the cloud form just a single sink 
particle. In this case an extended disk builds up around 
the star, since the seed turbulence gives rise to some low 
level of rotation in the collapsing core, but otherwise the 
gas contains little structure. 

Interestingly, there is no clear trend linking the number 
of fragments that form and the initial turbulent energy: 
the 0.4 c s run fragments more than the 0.8 c s run, de- 
spite containing only one quarter of the initial turbulent 
energy. The effects at play here are somewhat complex. 
First, the turbulent velocity field contained in the col- 
lapsing region will differ with each value of Au t urb/c s , 
since the different strengths of flow will push the gas 
around to different degrees. In addition, the nature of 
the turbulence that survives in the collapsing core will 
also affect the fragmentation. As we see from the im- 
ages in Fig. [T] the cloud with A^turb = 0.4c s - the most 
successful in terms of fragmentation - forms a large disk- 
like structure. As we will discuss in ^5] this configura- 
tion appears to aid the fragmentation of the infalling 
envelope. Therefore, much of the cloud's ability to frag- 
ment depends on the level of rotation that happens to 
become locked-up in the collapsing region. Further, and 
to a lesser degree, the extra delay in the collapse caused 
by the increased turbulent support also gives the cloud 
more time to wash out anisotropies in the gas (see Fig. [8] 
for the collapse times). Thus the ability of the cloud to 



fragment is a competition between these conflicting pro- 
cesses. For the randomly generated velocity field that we 
used in this study, the ability to fragment is better for 
Awturb/cs = 0.4, than Aw t urb/c s = 0.8, but we note that 
this may not always be the case. We stress that to make 
a quantitative statistical statement about the number of 
fragments that form as a function of the turbulent Mach 
number, we would need to run a series of different real- 
izations of the turbulent velocity field in each case. Such 
a comparison lies outside the scope of our current study. 

The mass functions of the sink particles from the sim- 
ulations that undergo fragmentation are shown in Fig. [2j 
For clarity, we have omitted the single 100 M sink parti- 
cle that forms in the 0.1 c s cloud. For the 0.2 c s and 0.8 c s 
clouds, we see that the sink masses cluster around some 
central value - roughly 12 M Q and 4 M Q respectively 
- while the 0.4 c s cloud has a mass function that is 
skewed to lower masses, with a peak at around 1 M Q , 
a sharp fall-off below this, and a broad distribution to- 
wards higher masses, extending up to around 13 M Q . 
We emphasize that the mass functions presented here (as 
well as in Fig. 10) do not yet represent the final IMFs, 
as they correspond to an intermediate time in the overall 
accretion process, where only 10% of the cloud has been 
accreted. They also might be affected by the numerical 
details of our sink technique, in particular the absence of 
any sink-sink mergers, as discussed in Section 3. 

The reason for the spread of masses becomes apparent 
when we look at the accretion properties of the sink par- 
ticles, which are shown in Fig. [3} Focusing on the 0.4 and 
0.8 c s runs, and looking at the evolution of the individual 
sinks in more detail, we sec that some appear to accrete 
rapidly and then suddenly stop. This behaviour is typ- 
ical of what is seen in simulations of bound, fragment- 
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Fig. 4. — Column density images showing the state of the 150 Mq Pop. III. 2 clouds after they have converted 10 percent of their mass 
(15 M©) into sinks. Note that the scale in this figure differs from that in Figfl] The clouds exhibit a different behaviour from their Pop. 
III.l counterparts. Although the clouds all form disks around their sink particles, due to the angular momentum contained in the initial 
turbulent motions, only the A«turb = 0.8c s turbulent cloud has undergone fragmentation by this point in the simulation. Note that the 
sink particles in this study have an accretion radius of 20 au. 



ing cores in the context of present-day star formation 



(Klessen k Burkert 2000 Klessen|2001 


Klessen k Burk- 


eTtl2001i Bonnell, Vine k Bate 20041137 


lmcja & Klessen 


20041, and is a result of velocity kicks from dynamical 



rate in such a system ( |Bonnell et al.pOOla ) is given by 

(6) 



m 

oc 

v 



where m* is the mass of the sink particle, p is the gas den- 
sity, and i> rc i is the velocity of the sink particle relative to 
the gas. The ability of a sink to accrete more mass from 
the available reservoir is significantly reduced once its 
velocity increases. The effect is then exacerbated by the 
fact that an increase in velocity results in the sink particle 
moving to a more distant orbit (or even being kicked out 
of the system entirely) and hence into a position where 



the gas density is lower. Since this sink particle is now in 
a position where further accretion is difficult, its siblings 
are able to accrete its 'share' of the mass reservoir, with 
the majority going to those few sinks that sit right in 
the middle of the cloud's potential well. In general, as 
the sinks accrete from the background gas, they tend to 
move towards the centre, due to mass-loading. Further, 
their increased mass makes them more likely to survive 
dynamical encounters with their less massive siblings. As 
such, the 'rich get richer', with a few sinks ending up sig- 
nificantly more massive than the rest. The process is 
typically termed 'competitive accret ion', (Bonnell et al. 



19971 |2001a|b} [Bonnell k Bate|2006| and normally leads 
to the type of distribution of masses seen in our 0.4c s run 



i type 
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Looking at the mass evolution of the individual sink 
particles, we also see that the formation of new sink par- 
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Fig. 5. — Temperature as a function of number density for the 
Pop. III.l (dark blue) and Pop. III.2 (light blue) Av tulb = 0.1 c s 
simulations. In both cases, the curves denote the state of the cloud 
at the point just before the formation of the sink particle. 



tides occurs in bursts. The turbulence generates struc- 
ture in the gas which is enhanced by the gravitational 
collapse as the gaseous envelope falls in towards the cen- 
tral system of sink particles. The bursts in sink formation 
reflect the moments when these structures detach from 
the flow and become self-gravitating in their own right. 

Figure [3] also shows the accretion rate of the cluster 
as a whole, and how that compares to an estimate made 
from the radial infall profile. To construct this estimate, 
we first compute the mass infall rate as a function of the 
radial distance r from the densest SPH particle: 



m(r) = Air r 2 p(r) v r (r), 



(7) 



where p(r) is the gas density in a spherical shell with ra- 
dius r and width dr, and v r (r) is the radial velocity of the 
gas in this shell and all quantities are volume-averaged. 
Given the enclosed mass as a function of radial distance, 
m enc (r), it is straightforward to convert from m(r) to an 
infall rate as a function of enclosed mass, m(m cnc ), which 
we identify with the total mass accretion rate of the sys- 
tem of sink particles. Figure [3] demonstrates that this 
estimate provides a fairly close match to the actual ac- 
cretion rate onto the sink particle population, except at 
very early times, which is a numerical artifact: the sink 
particles instantly accrete all gas within their accretion 
radius when they form. At later times, the main differ- 
ence between the estimated accretion rate and the true 
accretion rate is that the latter is significantly noisier, 
owing to the bursts of sink formation and the dumpiness 
of the infalling gas. Finally, we note that the accretion 
rates are also fairly insensitive to the level of turbulence 
in the cloud, suggesting that random turbulent motions 
with the magnitudes considered here should not lead to a 
variation in the overall accretion rates of primordial stars 
or star clusters from minihalo to minihalo. Other sources 
of support against gravity, such as rotation and magnetic 
fie lds (assuming the latter can be efficiently generate d, as 
in |Tan fc Blackman|2004| and |Schleicher et al.|20lo| , are 



likely to play a greater role in regulating the accretion 
rate. Indeed, the important role played by rotation in 
regulating the infall and accretion of gas can be appreci- 
ated if we compare the accretion rates measured in our 
simulations with those estimated or measured in previous 
studies of Pop. Ill star formation starting from more real- 



istic cosmological initial conditions (e.g. Abel, Bryan, & 
Norman|2002[|Bromm fe Loeb|2004| . The absence of ini- 
tial rotational support m our simulations leads to higher 
infall velocities, and hence to an accretion rate that is a 
factor of a few larger than these previous values. 

4.2. Pop. III. 2 clouds 
4.2.1. Small clouds with 150 M Q 

The first clouds we will examine in the Pop. III. 2 case 
are those with an initial energy balance similar to those 
studied in the Pop. III.l case: that is, clouds that have 
only a few Jeans masses initially. Since the Pop. III. 2 
channel is cooler at number densities around 10 5 cm~ 3 , 
with a typical temperature of around 75K, the same ini- 
tial number of Jeans masses requires that the clouds have 
a lower mass of 150 M©. 

The column density distribution in these clouds after 
10 percent of the gas has been accreted by the sink par- 
ticles is shown in the column density images in Fig. [4j 
We see from the images that these clouds undergo signif- 
icantly less fragmentation than their Pop. III.l counter- 
parts, forming at most 3 sink particles after 10 percent 
of the cloud mass has been accreted. While these calcu- 
lations demonstrate that this primordial star formation 
channel is susceptible to fragmentation if turbulence is 
present in the collapsing gas, it appears to be signifi- 
cantly more stable than the Pop. III.l channel that gives 
rise to the first stars in the universe. 

The main reason why this mode of star formation is less 
susceptible to fragmentation has to do with the thermal 
evolution of the gas as it collapses. In Fig. [5] we show 
the temperature of the g function of its number 

density, for the Pop. III.l and Pop. III. 2 cases, taken 
from the simulations with Au tur b — 0.1 c s . The elevated 
H2 and HD fractions in the Pop. III. 2 simulations do 
not provide enough cooling to keep the gas close to the 
CMB temperature at these densities, owing to the in- 
creasing inefficiency of the HD as a coolant as it nears 
the critical density at which the populations of its rota- 
tional level reach their local thermodynamic equilibrium 
(LTE) values. The gas therefore heats up significantly 
as it collapses, increasing its temperature from 75 K at 
n = 10 5 cm~ 3 to roughly 700 K at n = 10 7 cm -3 , cor- 
responding to evolution that is almost adiabatic. This 
sharp increase in temperature temporarily increases the 
local Jeans mass at the center of the cloud, and signif- 
icantly slows the collapse, giving the turbulence in the 
cloud time to decay. In the absence of additional physi- 
cal processes able to replenish the turbulence, there is 
nothing to sustain the density inhomogeneities in the 
cloud that act as the seeds for later fragmentation (since 
they are not yet self-gravitating), and so the end result 
is a collapse with a much lower level of fragmentation 
than in most of the Pop. III.l simulati ons. A similar 
effect has previously been noted by |Tsuribe fc Omukai] 
( 2008[ ) in their study of fragmentation in very metal poor 
gas clouds. They find that for metallicities of around 
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Fig. 6. — As Fig. [3] but for the 150 Mq Pop. III. 2 simulations. Again, the evolution is plotted until slightly more than 10 percent of the 
cloud's mass has been accreted, which in these cases occurs after roughly 600 yr. 
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Fig. 7. — Accretion rates as a function of enclosed gas mass 
in the Pop. III. 1 (upper lines; blue) and Pop. III. 2 (l ower lines; 
magenta) simulations, estimated as described in Section |4.1| Note 
that the sharp decline in the accretion rates for enclosed masses 
close to the initial cloud mass is an artifact of our problem setup; 
we would not expect to see this in a realistic Pop. Ill halo. 



Z ~ 10~ 45 Z Q , heat input due to three-body H 2 forma- 
tion at n ~ 10 8 cm~ 3 leads to a sharp jump in the gas 
temperature, which delays the collapse, reduces the elon- 
gation of the collapsing core, and suppresses any frag- 
menta tion. Furthermore, Yos hida, Omukai &: Hernquist| 
(2007) also briefly addressed this issue in their study of 
Pop. III. 2 star formation, and showed that their simu- 
lated Pop. III. 2 prestellar core would be stable against 
gravitational deformation at similar densities, owing to 
its hard effective equation of state. 
Although the evolution of the temperature with den- 
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Fig. 8. — Time taken for the first sink particle to form in each 
simulation. For reference, the horizontal dashed line denotes the 
free-fall time of the clouds at their initial density. 



sity is significantly different between the Pop. III.l and 
Pop. III. 2 channels, we see that the accretion rates are 
similar when we consider the inner 0.01 to 10 solar masses 
of the colla psing envelope (F ig. |6|, consistent w ith the 
results from Yoshida, Omukai & Hernquist ( 2007 1 . How- 
ever, if we consider the evolution of the accretion rates 
over the whole cloud, we see that at later times the two 



channels depart significantly from one another, with ac- 
cretion occurring at a significantly slower rate in the Pop. 
III. 2 case (Fig. |7]). Also, at very early times, the accre- 
tion rates in the Pop. III. 2 clouds are more sensitive to 
the level of turbulence than the Pop. III.l clouds, and the 
turbulence has a somewhat stronger effect in delaying the 
onset of star formation in these calculations (Fig. [M. 
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Fig. 9. — As Figs. ^and[4] but for the 1000 Mq Pop. III. 2 clouds after the sinks have accreted 10 percent of the total cloud mass (100 
Mq). Note that the scale in this figure differs from that in Fig. ^ Although these simulations have significantly more mass than those 
shown in Fig. [4] and are therefore initially more Jeans unstable, only a small amount of fragmentation is seen. The number of fragments 
formed in each case is much smaller than in the corresponding Pop. III. 1 runs. Note that the sink particles in this study have an accretion 
radius of 20 au. 



4.2.2. Large clouds with WOO M 

Given that the Pop. III. 2 clouds seem to be much more 
stable against fragmentation, it is worthwhile investigat- 
ing whether they can be made to fragment when the 
gas is initially more Jeans unstable. Recent simulations 
of the formation of the first galaxies show that regions 
where Pop. III. 2 star formation occurs are fed by cold 



supersonic turbule nt streams of gas (e.g. Greif et al. 2008 



Wise fe Abel|2008 ). As such, the initial condition for the 



Pop. 111.2 channel in this picture may have significantly 
more than one Jeans mass, due to the rapid assembly 
of the self-gravitating core. In this section we consider 
Pop. III. 2 clouds that contain 1000 M , but otherwise 
have the same properties as the 150 M clouds (i.e., same 
initial temperature, density, and initial chemical compo- 
sition). Since these clouds are colder than the Pop. III.l 
clouds of the same mass and density, they are initially 
more Jeans unstable, having around 26 Jeans masses at 
the start of the simulation. 

The results of two such calculations are shown in Fig.[9j 
in which the initial levels of turbulence have been set to 
0.4 c s and 0.8 c s . Again we show the evolution at the 
point where 10 percent of the cloud's mass has been con- 
verted into (or accreted onto) sink particles. Although 
these clouds are now initially more unstable than the 
Pop. III.l clouds, they form significantly fewer sink par- 
ticles, only 8 in the 0.4 c s cloud and 7 in the 0.8 c s cloud 



(see Fig. 10). The fact that the mass functions of the 



sink particles appear to have two peaks is due to the 
cloud undergoing two distinct bursts of sink formation. 

The interesting result here is that these clouds exhibit 
less fragmentation than their Pop. III.l counterparts, de- 
spite containing more Jeans masses in the initial config- 
uration. The main reason for this is the gas in the Pop. 
III. 2 clouds follows an extremely 'stiff' effective equa- 
tion of state (essentially adiabatic evolution), as can be 



seen from the temperature-density relationship shown in 
Fig. [5] Such a rapid increase in the temperature dur- 
ing the initial collapse makes generating structure in the 
gas difficult, helping to remove any anisotropies intro- 
duced by the turbulent flows. An additional effect is 
that for equations of state with effective adiabatic index 
of 7 > 4/3 (as is the case with these clouds), the Jeans 
mass increases with increasing density. As such, the col- 
lapse halts until sufficient mass has been assembled, and 
the new Jeans mass has been reached. Since the cooling 
time is significantly longer than the free-fall time, the 
gas has a chance to remove structure that could poten- 
tially assist the fragmentation at higher densities, where 
the effective 7 is more similar to that found in the Pop. 
III.l clouds. The combination of these effects results in 
a gas which is much more stable to fragmentation during 
collapse than in the Pop. III.l simulations. 

5. LONG TERM EVOLUTION OF THE INFALLING 
ENVELOPE 

More insight into the stability of the clouds against 
fragmentation can be gained from Fig. |11[ where we show 
the temperature and H2 fraction as a function of density, 
for the 1000 M Pop. III.l and III.2 clouds in which the 
initial level of turbulence was 0.4 c s . The different col- 
ors correspond to different points in the evolution of the 
clouds, with red, green and blue corresponding to the 
formation of the first sink particle, 50 M accreted, and 
100 M accreted, respectively. The temper ature-density 
plots show the sa me behavior as reported in Stacy, Greif] 
& Bromm (2010). The particle evolution appears to di- 
verge at a density of around 10 10 cm -3 , with one group of 
SPH particles heating up to a maximum temperature of 
around 7000K and the other staying around 1500K and 
cooling significantly at higher densities. As discussed in 
Stacy, Greif & Bromm (2010), the hot part of the di- 
agram corresponds to gas that falls in at a later stage 
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Fig. 10.— (top) Mass functions for our two 1000 M Q Pop. III. 2 
clouds. As with the Pop. III. 1 data in Fig. [2] we show the mass 
function after 100 Mq of gas has been accretedoy the sink particles. 
Due to the reduced fragmentation in these calculations, the sink 
particles are on average more massive than their Pop. III. 1 coun- 
terparts since they have less competition for the available mass, 
(bottom) The combined mass function for the two simulations. 



in the evolution of the cloud. At that point, the en- 
closed gas mass is larger than at earlier times, and so 
the free-fall velocity is correspondingly larger. The gas 
therefore shocks more strongly than at early times, caus- 
ing it to become hot enough to collisionally dissociate its 
H 2 rapidly. With the H 2 gone, there is nothing to cool 
the gas until its temperature reaches T ~ 7000 K. At 
this point, Lyman-a cooling becomes effective, allowing 
it to resist further heating. Although we see the same 
trends in both the Pop. III. 1 and Pop. III. 2 cases, in the 
latter the amount of gas departing from the 'standard' 
temperature-density evolution is significantly reduced, 
since the elevated H 2 fractions allow the post-shock gas 
to cool more effectively limiting the temperature rise and 
allowing more of the gas to retain its H2. 

From the accompanying graphs we see that the H 2 is 
rapidly dissociated as the temperature rises, but given 
the relatively low amounts of H 2 present in the gas at 
these densities, the cooling provided by the dissociation 
is clearly unable to offset the compression. However we 
see that at higher densities the H 2 suddenly reforms, co- 
inciding with the region in the temperature-density dia- 



grams where the gas is cold. This prompts the question 
of whether the cold gas is the result of an isochoric cool- 
ing instability, brought on by rapid H 2 formation? Does 
the turbulence trigger an instability that naturally exists 
in pure primordial gas, but that is less effective in gas 
that has been influenced by the presence of previous star 
formation? 

By harnessing the Lagrangian nature of SPH, the plots 
in Fig. [12] shed some light on this issue. They show the 
temporal evolution of the density, temperature and H 2 
fraction for several SPH particles that eventually become 
sink particles, and as such trace the conditions in the 
gas in the run-up to gravitational fragmentation and col- 
lapse. Again the simulations are those shown in Fig. 11 



The quantities are calculated by averaging over each par- 
ticle's 50 nearest neighbors at each instant in time. In 
the case of the Pop. III. 1 cloud, these sinks are the last to 
form in the simulation, while in the Pop. III. 2 cloud, the 
five lines represent all of the sink formation that occurs 
between the formation of the first and last sinks. 

The figure shows a number of interesting features. 
First, we see that none of those SPH particles destined to 
become sinks undergo the rapid rise in temperature - and 
accompanying loss of H 2 - that is shown in Fig. [TT] In 
contrast, their temperatures remain close to the temper- 
atures found within the first collapsing core (see Fig. IB"]). 
This demonstrates that the cool particles seen in Fig.]ll| 
do not come from regions that undergo shock heating and 
subsequent loss of H 2 . Instead, we see that they come 
from regions of gas that first undergo a relatively quies- 
cent collapse, before being involved in several expansions 
and contractions. 

The fact that this first stage of the density evolution 
is fairly slow needs to be stressed: the free-fall times at 
densities of n = 10 9 cm~ 3 and n = 10 10 cm -3 are ap- 
proximately 1600 yr and 500 yr, respectively. As such, 
these particles are not experiencing as much compres- 
sion as those that end up losing their H 2 . In fact we see 
that they actually have a higher than average H 2 frac- 
tion, when we compare them with the lower right-hand 
plot in Fig. |11[ a property that helps them remain fairly 
cool as they collapse. The reason why they experience 
less compression is that they are collapsing in a rotating 
structure that has been formed during the collapse of the 
turbulent gas. 

The evolution of these particles and their immediate 
surroundings alters abruptly once they enter the dynam- 
ically complicated swirling regions that we can see in the 
column density figures. First, they collide with other ma- 
terial, which results in a sharp increase in their temper- 
ature and density. This increase in the density in turn 
increases the rate of H 2 formation, which very rapidly 
turns them fully molecular. As they re-expand, the adia- 
batic cooling and the now significantly enhanced H 2 line- 
cooling act together to reduce the temperature. In some 
cases this happens several times, but for all of these par- 
ticles, the end result is the same and they find themselves 
in a clump of gas that is now Jeans unstable and fully 
molecular: the perfect conditions for forming a new pro- 
tostar. Interestingly, we see similar behaviour leading up 
to the formation of the sink particles in both the Pop. 
III. 1 and III. 2 clouds. The relative lack of fragmenta- 
tion in the Pop. III. 2 case simply results from the lack 
of structure in the collapsing envelope, rather than any 



13 




n [cm 3 ] n [cm 3 ] 




n [cm 3 ] n [cm 3 ] 

Fig. 11. — Temperature and H2 fraction as a function of density in the Pop. III.l (top) and 1000 Mq Pop. III. 2 (bottom) simulations in 
which Ai> tur b = 0.4c s . Three stages in the evolution of the clouds are shown, corresponding to when the first sink forms (red), when there 
is 50 Mq in sink particles (green) and finally when there is 100 Mq in sink particles (blue). 



thermal properties of the gas at these high densities. 
6. DISCUSSION 

The calculations presented in this study suggest that 
the first stars in the universe may form in small dense 
clusters, provided that the turbulent initial conditions 
we adopt are close to those found in minihalos. This 
suggests that the Population III IMF covered a broad 
range in masses, possibly exhibiting a scale-free, power- 
law extension similar to the present-day case. Previous 
arguments in favor of a peaked IMF, in the shape of a 
narrow Gaussian or even a delta function, would then 
need to be revisited (e.g., Bromm & Larson 2004). Al- 
though a similar prediction can be drawn from the study 
of Clark et al. (2008, hereafter CGK08), there are some 
subtle differences. Our present study employs a fully 
self-consistent treatment of the thermodynamics, rather 
than the piece-wise polytropic equation of state (EOS) 
approach used by CGK08, which was taken as a fi t to th e 
detailed one-zone calculations of jOmukai et al.| (20051. 
Although both start with similar (Pop. 111.1) initial con- 
ditions, the fragmentation seen in the current simulations 
has a different origin. In the self- consistent treatment 
employed here, the fragmentation is driven by the com- 
plicated thermodynamics of high density clumpy gas as 
it enters the disk-like regions that surround the first pro- 
tostar. In contrast, the fragmentation in CGK08 was 



evident at much lower densities, with the turbulent flows 
forming structures that were enhanced during the col- 
lapse, and rotation providing a 'window of opportunity' 
for those structures to become gravitationally unstable in 
their own right, rather than being simply accreted onto 
the central protostellar core. Note however that the sim- 
ulations in CGK08 also contained systematic, solid-body 
rotation in the initial conditions, which we do not study 
in this paper. This may help some of the structure to 
survive. 

In general, the calculations presented in this paper sug- 
gest that the structure formed early on in the collapse is 
less likely to survive when the full thermodynamical be- 
havior of the gas is taken into account. This is evident in 
the fact that the current calculations yield significantly 
fewer fragments than the Pop III simulations in CGK08 
- 25 sink particles for 19 M Q of accreted gas in CGK08, 
compared to 31 sink particles for 100 M of accreted gas 
in the current study. One feature that sh ould be stressed 
(and which w as originally pointed out in Stacy, Greif & 
Bromm|2010 ) is that the long-term thermodynamic evo 



lution of the envelope differs significantly from the one- 
zone models. This brings into question the practice of 
using a piece- wise barotropic EOS as a proxy for the full 
thermodynamics in simulations that intend to study the 
evolution of the gas beyond the collapse of the first pro- 
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Fig. 12. — Particle trajectories showing the physical conditions in the gas leading up to the formation of sink particles in the 1000 
Mq Pop. III. 1 and Pop. III. 2 clouds, from the simulations with A« tur b = 0.4 c s . The quantities are calculated by averaging the density, 
temperature and H2 over the 50 nearest neighbours of the particles which eventually become sink particles. Their evolution is shown from 
the onset of the star formation in the cloud (i.e. from the formation of the first sink particle), up to the point at which they themselves are 
turned into sink particles. 



tostellar core. 

Assuming that the turbulent initial conditions used in 
this paper are indeed representative of the gas in miniha- 
los, it is worth considering the implications of the frag- 
mentation that we see. The mass spectrum of the frag- 
ments in our calculations ranges from a few 0.1 Mq to 
a few 10 Mq. We would not expect stellar feedback 
to chan ge this result signific antly, since previous stud- 
ies (e.g. McKee fc Tan|[2008| ) have shown that feedback 
becomes effective at limiting accretion onto Pop. Ill pro- 
tostars only for protostellar masses greater than about 
20 M . However, such feedback effects are expected to 
become important, and in particular to suppress further 
fragmentation (e.g., Krumholz et al. 2009) during the 
further evolution of the cloud. The fragments will then 
likely grow in mass, and possibly even merge with each 
other. It is very difficult to extrapolate to the final situ- 
ation where accretion and merging stops. But our main 
result that the Population III IMF was likely broad seems 
robust. 

As we have discussed above, it is difficult to reach any 
definitive conclusions regarding the final Population III 
masses and the resulting IMF. It appears at least possi- 
ble, however, that in rare cases truly metal-free stars with 
masses less than ^0.8 M Q could have formed, and would 
still be present in the Milky Way today, thus providing a 
unique opportunity to directly probe the physical condi- 
tions at the end of the dark ages with Galactic observa- 
tions. Current models of hierarchical galaxy formation 
predict that the first and most metal-poor halos to merge 
will become part of t he bulge compone nt of the result- 
ing spiral galaxy (e.g. |Tumlinson]|2010[ ) . It is therefore 
highly interesting to survey the bulge of the Milky Way 



for extremely metal-poor and metal-free stars. However, 
this is also very challenging as the bulge is far away, has 
high stellar density, and contains a wide range of stellar 
populations of all ages and metallicities. The odds of 
finding a few truly metal-free stars amongst millions of 
other stars are very low, or even zero due to pollution 
from the interstellar medium (Frebel et al. 2009). 

With the same caveat as above, it is interesting to con- 
sider the implications of a Pop. Ill IMF with a char- 
acteristic mass significantly smaller than the canonical 
100 Mq. For instance, the flux of ionizing photons pro- 
duced by a population of Pop. Ill stars has a significant 
dependence on the form of the Pop. Ill IMF, as does 
the metal-enrichment p attern produced by a collection 
of Po p. Ill supernovae (Tumlinson, Venkatesan & Shull 
2004|. Indeed, the abundance patterns observed to date 



m extremely metal-de ficient stars in the Galac tic halo 
(see e.g. the review by Beers & Christlieb 2005) are far 
more consistent with an IMF that produces primarily 
core-collapse supernovae, with progenitor masses of 10- 
40 Mq , rather than with an IMF that produces only very 
massive pa ir-instability supernovae, or PISNe (Joggerst 
et al.||2010 1 . On the other hand, there may be subtle se 



lection effects at work that bias current surveys against 
finding PISN-enriched stars (Karlsson et al. 2008). The 
basic argument here is that PISNe have such high metal 
yields that abundances in stars that form out of this ma- 
terial are already quite high (see Greif et al. 2010), and 
would therefore be missed in searches that target the low- 
est metallicities. The interpretation of the large carbon 
enhancements seen in the population of carbon-enhanced 
metal-poor stars as the result of the enrichment of these 
stars by winds from binary companions that have passed 
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through the AGB phase also implies the existence of a 
large number of inter mediate-mass Pop. HI stars, with 
masses M = 1-8 M ( |Tumlinsori]|2007a|b[ ) . Again, there 
are alternative models to explain the carbon enhance- 
ment in metal-poor stars in terms of nucleosynthesis in 
faint supernovae (e.g., Iwamoto et al. 2005), in line with 
a higher characteristic Pop. Ill mass. 

It is important to note that one factor that may limit 
the impact of turbulent fragmentation on the primor- 
dial IMF is if most Pop. Ill stars form in conditions re- 
sembling our Pop. III. 2 clouds. It is relatively straight- 
forward to show that most Pop. Ill stars will form in 
halos that have been affected in some f ashion by a previ- 
ous episode of Pop. Ill star formation dTrenti k. Stiavelli 
~ Greif et al.J |2010 |>. Thus, using the |Tan fc McKee 



2009 



I 2008) terminology, most Pop. Ill stars will be Pop. 111.2 



stars. However, it is less obvious how many of these stars 
will have formed out of gas that has been cooled to tem- 
peratures T -C 200 K by HD. The effectiveness of HD 
cooling in these systems depends on the balance between 
the enhanced formation of H2 and HD, owing to the en- 
hanced initial fractional ionization in this gas, and the 
destruction of H 2 and HD due to Lyman- Werner band 
absorption of ultraviolet photons from an extragalactic 
background (Haiman, Abel & Rees 200 0; Johnson, Greif 
& Br omm 2008), or from lo cal sources ( |Omukai 



Nishi 



1999 Glover fc Brand|2001 ). The relative importance of 



these effects has not been explored in great detail, and in 
any case, the outcome is likely to be sensitive to uncer 
taint ies i n the microphysics that have on l y recently been 



resolved (Glover, Savin fc Jappse n 2006 Glover & Abel 
20081 |Kreckel et al.||20lo[ )~ 

To sum up, our calculations suggest that turbulent 
fragmentation may play an important role in the forma- 
tion of Pop. Ill stars, and may strongly influence the form 
of the Population III IMF. However, the approximate na- 
ture of our calculation - specifically, our simplified choice 
of initial conditions - means that we should regard this at 
present as no more than a plausible hypothesis. To bet- 
ter establish the role of turbulent fragmentation in real 
Pop. Ill minihalos, we will need to be able to simulate a 
representative sample of both Pop. III. 1 and Pop. III. 2 
minihalos with sufficiently high resolution, such that they 
can resolve the turbulent flows in the gas at the moment 
when it becomes self-gravitating. This is a challenging 
prospect, lying far beyond of the scope of this prelimi- 
nary study, but our results suggest that it would prove 
extremely worthwhile. 

7. SUMMARY 

We have explored the effects of subsonic turbulence on 
the gravitational collapse of primordial gas clouds. The 
study employed sink particles to model the run-away col- 
lapse of protostellar cores, which allowed us to follow the 
evolution of the collapsing clouds beyond the formation 
of the first protostar. The calculations also used a full 
time-dependent chemical network that accounts for the 
thermodynamic behaviour of the gas. The current calcu- 
lations contain neither magnetic fields nor feedback from 
the protostars. Our main findings can be summarized as 
follows: 

• Turbulent primordial gas is unstable to fragmenta- 
tion when one considers the evolution beyond the 



formation of the first protostellar core. 

• Gas starting from conditions appropriate to Pop. 
III. 1 collapse, rather than Pop. III. 2 star forma- 
tion, is m ore susceptible to fragmentation. As sug- 
gested by Yoshida, Omukai & Hernquist ( 2007 ), the 
thermal evolution of Pop. 111.2 gas as it collapses 
helps to suppress further gravitational instability 
over and above the main collapse mode. However 
some fragmentation in the Pop. III. 2 case is seen 
in this study, caused by the inhomogeneities intro- 
duced by the turbulence. 

• In the cases where fragmentation is efficient (in par- 
ticular the Pop III. 1 cloud with turbulent rms ve- 
locity Aw turb = 0.4c s ), the masses of the fragments 
extend over a large range, which results in a dis- 
tribution of stellar masses exhibiting a power-law 
extension towards high mass, as seen in the present- 



day IMF ( [Kroupa |2002[ |Chabrier| 2003| . However 
it should be noted that the exact form of the sink 
particle mass function may depend on our assumed 
initial conditions (see Appendix [C]) and the simpli- 
fications employed in the sink particle implemen- 
tation itself. 

• Due to the relative lack of fragmentation in the 
Pop. III. 2 clouds compared to the Pop. III.l clouds, 
the stars formed in the Pop. III. 2 calculations 
are of higher mass on average (for the same ac- 
creted mass) than their Pop. III.l counterparts, 
even though the mass accretion rate of the star 
cluster as a whole is higher in the Pop. III.l clouds. 

• Fragmentation tends to occur in gas which has 
temperatures of around 200 - 400 K at densities 
above 10 11 cm -3 , significantly lower than the 1000 
- 1500 K associated with the collapse of the first 
core. Such cold temperatures are a result of the ex- 
pansion that occurs as gas enters the rotating, disk- 
like regions around the central core, coupled with 
relatively high fractions of H2, which can provide 
efficient line-cooling. These rotating structures are 
themselves a consequence of the angular momen- 
tum that is present in the turbulence. 

In summary, we propose that if even small levels of 
turbulence, with velocity dispersions of order 20% of the 
sound speed or more, are present in the baryonic compo- 
nent of dark matter minihalos, primordial stars are likely 
to be born in small stellar groups, rather than in isola- 
tion, and to have a wide range of stellar masses. Further, 
the very first stars (Pop. III.l) may have lower masses on 
average than the second generation of stars (Pop. III. 2), 
contrary to what has previously been assumed. 
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APPENDIX 



CHEMISTRY 

To model the chemical evolution of the metal-free gas, we use the network for primordial hydrogen, helium and 
deuterium che mistry detailed in Table |A"lj Our treatment of the hydrogen and helium chemistry largely follows [Glover| 
& Abel (2008), but our treatment of the deuterium chemistry is significantly simplified. This simplification arises from 
our neglect of D - , HD + and D2, none of which play a significant role in controlling the HD abundance in the physical 
condit ions relevant to our present study. For the most part, our choice of rate coefficients also follows [Glover fc Abel| 



(I2OOSJ) . 2 We adopt case B rat e coefficients fo r the r ecombination of H + and He ++ , and treat He + recom bination 
as decribed in Section 2.1.4 of Glover & Abel (20081. We adopt rate coefficients from Galli & Palla (1998) for the 
associative detachment of H~ ions by atomic hydrogen (reaction 2), and for the mutual neutralizatio n of H by H + 
(react ion 5). We also note that the uncertainties in these reaction rate coefficients discussed in |Glover, Savin fc Jappseri 
( |2006[ ) are unlikely to be significant in the high density, low ionization gas modelled in our present simulations. For 
the rate coefficient of reaction 30, the three-body formation of H2 with atomic hydrogen as the third body 

H + H + H^H 2 + H, (Al) 

we adopt the rate coefficient proposed b y 
uncertain (|Glover||2008| iTurk et al.|2010' 



Glover (2008). The value of the rate coefficient for this reaction is highly 
jut our chosen value lies intermediate between the fastest and slowest rates 



to be found in the literature ( |Flower fc Harris 2007| and Abel, Bryan, fc Norman||2002| respectively), and is based 
on the best currently available data for the inverse reaction (i.e. collisional dissociation of H2 by H ; react ion 9). The 
influence of the uncertainty in this rate coefficient has been studied in detail elsewhere ( Turk et al.||2010 ). To fix the 
rate coefficient for the analogous reaction 

H + H + H 2 ^H 2 + H 2 , (A2) 



we follow Palla, Salpeter, fc Stahler (19831) and assume that the rate coefficient is one-eighth of the size of the rate 
coefficient adopted to r reaction |A1| For Three-body formation in whic h neutral helium is the third b ody, we follow 



Glover fc Abe l (2008) and use a rate coefficient originally taken from iWalkauskas & Kaufmanl (|1975|) . We assume, 
following Flower & Harris ( 2007[ ), that the rate coefficients for the three- body formation of HD (reactions 43-45) are 
the same as those for the corresponding H 2 formation reactions (nos. 30-32). 

We adopt values for the rate coefficients for collisional dissociation of H 2 by H, H 2 and He that are consistent with 
our choices for the three-body association rates in the sense that each pair of rate coefficients individually satisfies 



^form 
^dest 



(A3) 



where k[ orm is the rate coefficient for three-body association, fcdes t is the rate coefficient for collisional dissociation, and 
where the equilibrium constant K is given in all three cases by ( |Flower fc Ha^[2007l ) 

/ 52000 \ 



K 



1.05 x 10" 22 T 



-0.515 



exp 



V T 



(A4) 



This procedure is necessary in order to ensure the correct chemical behaviour of the gas at high densities and tem- 
peratures, where the H 2 formation and destruction timescales are both short, and most of the gas is close to chemical 
equilibrium. 

N ote th at there are two typographical errors that we are aware of in the set of reaction rate coefficients listed in Table Al in |Glover fc| 



Abel 



(2008|l. First, the fitting function used to describe the rate coefficient for the reaction H2 +H+ — ¥ Hj +H (reaction number 7 in their 



table) should use natural logarithms, and not base 10 logarithms, as listed. However, the listed fitting coefficients are correct. Second, the 
temperature dependence of the reaction H2 + He — > He + H + H + (their reaction 24) should be exp(— 35/T), and not exp(+35/T). 
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TABLE Al 

Reactions included in our chemical model 



No. 




Reaction 


References 


1 


H + 


c~ — > H~ + 7 


1 


2 


H- 


4-H -> H 2 +e" 


2 


3 


H + 


H+ -> +7 


3 


4 


H + 


H+ -4 H 2 + H+ 


1 


5 


H- 


4- H+ -> H + H 


2 


6 


H 2 + 


4-e- -» H + H 


5 


7 


H 2 - 


h H+ -> H+ + H 


6 


8 


H 2 - 


r c - -> H + H + e~ 


7 


9 


H 2 - 


rH->H + H + H 


8 


10 


H 2 - 


- H 2 — > H 2 + H + H 


9, 10 


11 


H 2 - 


- He -» H + H + He 


11 


12 


H + 


c~ -> H+ + e - + e - 


12 


13 


H+ 


4-e - -> H + 7 


13 


14 


H- 


+■ e - ->H + e"+e" 


12 


15 


H~ 


+ H^H4H4e" 


12 


16 


H+ 


4 H" -5- H+ 4- e- 

- e~ — > He+ 4- e~ 4- e~ 


14 


17 


He 4 


12 


18 


He+ 


4- e~ -> He++ 4- e - 4- e - 


12 


19 


He+ 


4- e~ — > He 4- 7 


15, 16 


20 


He++ + e - -> He+ + 7 


13 


21 


H~ 


4- H+ -> H 2 + H 


17 


22 


H~ 


4HJ->H4H4H 


17 


23 


H 2 - 


he" 4H"4H 


18 


24 


H 2 - 


- He+ — > He 4- H 4- H+ 


19 


25 


H 2 - 


- He+ -s. H+ 4- He 


19 


26 


He+ 


4- H -4 He 4- H+ 


20 




He 4 


- H+ -5- He+ 4- H 


21 




He+ 


4- H- -> He 4- H 


99 


9Q 


He 4 


- H- -> He 4- H 4- e" 


93 


30 


H + 


H 4 H -> H2 4 H 


24 


31 


H + 


H 4- H 2 -4 H 2 4- H 2 


25 


32 


H + 


H 4- He -> H 2 + He 


26 


33 


D+ 


4-e - -> D 4- 7 


27 


34 


D + 


H+ -> H4-D+ 


28 


35 


H + 


D+ -> D4-H+ 


28 


36 


H 2 - 


r D^HD4H 


29 


37 


H 2 - 


-D+->HD4H+ 


30 


38 


HD 


4- H -> H 2 + D 


31 


39 


HD 


4- H+ -> H 2 4- D+ 


30 


40 


D + 


e - -> D+ 4-e~ 4- e _ 


27 


41 


He+ 


4- D -4 D+ 4- He 


32 


42 


He 4 


-D+^D4 He+ 


32 


43 


D + 


H4H->HD4H 


Sec text 


11 


D + 


H4H2 -+ HD + H 2 


Sec text 


45 


D + 


H 4- He ^ HD 4- He 


Sec text 



Not e. — 1: [Wishartl jl979t, 2: IGalli fc PaUal 
(1998), 3: Ramaker & Peck (1 976), 4:^Karpas, AnP] 
cich & Huntress (1979). 5: Schneider ct al. (1994), 6: 
Savin ct al. (2004), 7: Trcvisan & Tennyson (2002a), 
8: [ Martin, \jchwarz fc~K landv (19 96t, 9: [Martm, | 
Keogh Hi Mandy (1998), 10: Shapiro & Rang (1987), 
11: |l)ove et al.'l p^STf , 1 2: [janev el al-H 19871 , 13: 
Fcrland ct al. (1992), 14: Poulacrt ct al. (1978), 15: 
Hummer & Storey (1998), 16: Aldrovandi & Pcquig- 
not Hl973j , 17:j Dalgarno fc Lepp|(|1987^ 18: E chulz 
& Asundi (1967), 19: Barlow (1984), 20: Zygclman 
et al.| (|l989t, 21 : |Kimura efall j|1993) , 22: IFeart 
& MaytonHl994||, 237 | Muq et al.] | |l982| , 24: |Glover 
( 2008T T25: | (ilover | pOU8|, rescaled as in lllqwer fc 
|Harris| |200"7| , 26: IWalkauskas fc Kaufman| jl975fr7 
27: Same as corresponding 1 1 reaction, 28 : IKavinl 
IgO OgJ, 29: Fit t o dat a from IMielke et a l] fl2UU3B, 
3DTTC r erlich1 <1982t, 31: |ShavittHl959l, 32: iGlover fcl 

lAbJjfnnsV — ' 1 11 — ' 1 ' 
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COOLING FUNCTION AND THERMODYNAMICS 



The dominant coolant in our simulations is molecular hyd rogen. To mode l rotat ional and vibrational line emission 
from H2, we use the detailed cooling function described in Glover & Abel (20081, that includes contributions from 



collisions of H2 with H, He, H2, protons and electrons. At densities n £ 10 a cm 3 , the strongest of the H2 lines become 
optically thick, reducing its effectiveness as a coolant. To model H 2 cooling in this regime, we use an approach b ased 



on the Sobolev approximation that was first used in models of primordial star formation by Yoshida et al. 
write the H 2 cooling rate in optically thick gas as 



(2006). We 



Ah 2 = A£; ul A ul /3osc,ui^u, 



(Bl) 



u,l 



where n u is the number density of hydrogen molecules in upper energy level u, A.E U \ is the energy difference between 
this upper level and a lower level I, A u \ is the spontaneous radiative transition rate for transitions between u and I, and 
/?csc.ui is the escape probability associated with this transition, i.e. the probability that the emitted photon can escape 
from the region of interest. We take values for the level ene rgies from the compilatio n made available by P. G. Martin 



on his website 3 and for the radiative transition rates from Wolniewicz et al. (19981, and we fix n u by assuming that 
the H2 level populations are in LTE. The problem of modellin g op tically thic k H2 c ooling thereby reduces to one of 
computing the escape probability for each transition. We follow Yoshida et al. (2006) and write the escape probability 
for the transition u — > 1 as 

/?esc,ul = it (B2) 



where we approximate r u i as 



t u1 = a u iL a (B3) 

where a u \ is the line absorption coefficient and L s is the Soloblev length. In the classical, one-dimensional spherically 
symmetric case, the Sobolev length is given by 

Vth 



|di> r /dr 



(B4) 



where Vth is the thermal velocity, and dv T /dr is the radial velocity gradient. In our inherently three-dimensional flow, 
we generalize this as 



following |Neufeld fc Kaufman| ( |1993|) 
much larger than the size of the collapsing core 



If the velocity dispersion of the gas is very small, then L s can become very large, 
To ensure that we do not reduce the H2 cooling rate in this case to 



an artificially low value, we take as our actual length scale in Equation B3 the smallest of the Soloblev length and the 
local Jeans length, Lj. 

Since the line absorption coefficient a u i is linearly proportional to the number density of H2, we can write r u i as 



Tul 



Qui 
riH 2 



A H , 



,cff , 



(B6) 



where -/VH 2 ,cff = nn 2 L a is an effective H2 column density, and where a u \/nn 2 is a function only of temperature. We 
therefore tabulate the cooling rate per H2 molecule in the optically thick limit as a function of two parameters: the 
gas temperature T and the effective H2 column density AH 2l eff> and compute cooling rates during the simulations by 
interpolation from a pre-generated look-up table. 

At densities n > 10 14 cm -3 , a second form of H 2 cooling becomes important, called collision-induced emission. 
Although H2 molecules have no electric dipole, the interacting pair in a collision of H2 with H, He or H2 briefly acts as 
a 'supermolecule' with an non-zero electric dipole, and a hence a non-zero probability of emitting or absorbing a photon 
through a dipole transition. Because the collision time is very short, the re sulting collision-in duced transition lines are 
very broad, effectively merging into a continuu m (for more details, see e.g. Frommhold 1993"). We include this process 
in our cooling function, using a rate taken from Ripamonti & Abel (2004 1. We crudely account for the reduction of the 
CIE cooling rate by continuum absorption at very high number densities using the following prescription (M. Turk, 
private communication) 



where 



1 — eT TCIE 

AciE, thick = AciE.thin x niin ( , 1 

T~CIE 



T"CIE = 



(B7) 



(B8) 



However, we note that optical depth effects do not strongly affect the CIE cooling rate for densities below our threshold 
for sink particle creation (see £^3j) , and so the approximate nature of this opacity cutoff is unlikely to significantly affect 
our results. 



3 http:/ /www. cita.utoronto.ca/^pgmartin/h2. html 
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In addition to H2 ro-vibrational line emission and CIE cooling, our cooling function contains a number of other 
radiative processes: electronic excitation of H, He and He + , cooling from the recombinati on of H + and He + , Compton 
cooling and bremsstrahlung. Details of our treatment of these processes can be found in Glover & Jappsen (2007). 

We also account for changes in the thermal energy of the gas due to changes in its chemical makeup. Specifically, 
we include the effects of cooling due to the collisional ionization of H, He and He + , and due to the destruction of H2 
by charge transfer and by collisional dissociation, as well as heating due to the three-body H 2 formation. The balance 
between cooling due to H2 collisional dissociation and heating due to three-body H2 formation plays a very important 
role in regulating the thermal evolution of the gas at densities n > 10 s cm~ 3 . 

EVOLUTION OF THE BONNOR-EBERT SPHERES 

Given our somewhat arbitrary choice of initial conditions, it is prudent to ask whether they are applicable to 
primordial star formation, and in particular, whether the choice of the initial density profile affects the details of the 
collapse. In Fig. C13| we show, for the Pop III. 1 clouds with 0.1 and 0.4 c s turbulence, the number density and 
enclosed mass as a function of radius and how the radial velocity and temperature varies as a function of the enclosed 
mass. In the case of the 0.1c s cloud, the radial profiles (mas s and number density ) are much like those presented for 
the standard Pop III studies ( Abel, Bryan, & Norman 2002 Yoshida et al. 2006 1 , in which fully cosmological initial 
conditions were used. In particular, by the time that the first sink particle forms in this run, the gas has developed the 
same n oc r~ 2 ' 2 radial density profile as found in the fully cosmological runs, demonstrating that the density profile at 
this stage is insensitive to the shape of the initial density profile. 

We do however see differences between the radial velocity and temperature profiles, compared to those published in 
the literature. The infall velocities at this stage are somewhat higher than those seen in fully cosmological studies, 
but it should be noted that we do not include systemic rotation in our study, and so a major source of support is 
missing. However, aside from the higher values of the radial velocity, the shape of the profile is also different. In the 
runs presented in this study, the infall velocities of the shells enclosing 1 M and greater are systematically higher 
than those seen in the cosmological simulations published to date, with a peak in the radial velocity profile at about 
10-100 M (depending on the level of turbulent support) instead of around 0.3 M . Some of this difference in the 
radial velocity profiles is reflected in the somewhat steeper than reported temperature profiles, since the gas will reacts 
to the increased pdV heating. However the temperature profiles are also affected by the differences in the chemistry, 
and in particular, our choice of the three-body H 2 formation rate. 

So the obvious question is, do our high infall velocities make fragmentation more likely than it would be reality? 
To address this, we performed another Pop III. 1 simulation with 0.4 c s turbulence, but this time we set the confining 
external pressure to a factor of 4 less than was used in the original simulations. The radial a nd en closed mass profiles 

(right-hand panel). 



C13 



at the point just before the first sink particle is formed are shown as the dashed lines in Fig 
We see that radial velocity profile is affected by the change in boundary pressure, resulting in a factor of 4 reduction 
in the infall speed. Compared to the original run, which formed 31 sink particles, this simulation forms 22 sinks when 
10 percent of the gas has been accreted. While this suggests that the fragmentation we see is sensitive to the details 
of the collapse, and we are likely overestimating the number of fragments formed (at least at the scale of our sink 
particles) when one compares to the standard Pop III collapse profiles, it seems that this is not the main cause of 
the fragmentation. The gas still fragments, and the processes leading up to the fragmentation are the same as those 
discussed in f|5] 

Finally, we note that there is also new evidence to suggest that the radial collapse profiles seen in the literature 
so far, may n ot be representative of all Pop I II st ar formation. For exa mple, while the radial density profiles from 
the studies by Abel, Bryan, fc Norman ( 2002|) a nd Yoshida et al. ( 2006) are similar to the newer, higher-resolution, 
cosmological simulations of Turk, Abel & O'Shea ( 2009[ ) and Turk et al. ] |2010 1 , the radial (or enclosed mass) profiles of 
the other quantities in the n ewer calculations show considerable variati on. In particular, the velocity and temperature 
profiles in |Turk et al. (20101 lie somewhere in between those found by Yoshida et al. (2006) and our 0.4 c s turbulent 
cloud. A very high-resolution study of the onset of the initial gravitational instability in the baryon ic component of 
dark ma t ter m iniha los, may help to establish w hether these differences between the studies of, say, |Abel, Bryan, fc 
Norman ( 2002[ ) and Turk, Abel fc O'Shea||2009| are due to the increased resolution and improved physical treatment 
used in the latter study (e.g. the inclusion of the effects of three-body H2 formation heating), or just reflect the effects 
of cosmic variance. 
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Fig. C13. — Radial profiles from the Pop III. 1 clouds with A«turb = 0.1c s (left panel) and A« tar i, = 0.4c s (right panel). The lines show 
the initial conditions (solid line), the state of the gas just before the creation of the first sink (dotted line). In the case of the A% ur b = 0.4c s 
cloud, we also show the results of a collapse that has a factor of 4 less external pressure (dashed l ine). The dot-dashed line shows a slope 
of n a r~ 2 2 . The crosses and circular points denote approximate values taken, respectively, from |Yoshida et el!.| ( |2006[ | and |Abel, Bryan^] 
|fc Norman| p002l > 



